//!  追赶法求解三对角阵 (矩阵元迭代)
# ifndef TDMA_H_
# define TDMA_H_

#include <iostream>
#include <cmath>
#include <complex>

#define complex std::complex<double>


class TDMA {
private:
  int n; //n =nr
  complex *cs, *ds;
  //complex *a0, *a1, *a2;
  /*coefficient matrix A: 
  Subdiagonal Elements(a0), 
  Main Diagonal Elements(a1), 
  Superdiagonal Elements(a2)*/
  //constant matrix elements(d)

public:  
  //AX=D,X=LU,LY=D,  UX=Y
 
  TDMA (int n_): n(n_) {
    /*Obtain the subdiagonal elements of the upper triangular matrix U 
    from the coefficient matrix A*/
    cs = new complex [n-1];//U 
    ds = new complex [n];  //Y  

  }
  virtual ~TDMA () {
    delete[] cs;
    delete[] ds;
 
  }
  //AX=D
  template<typename T0, typename T1, typename T2, typename T>
  void solve (T0 a0, T1 a1,T2 a2, T d, T X1);

};

template<typename T0, typename T1, typename T2, typename T>
void TDMA::solve (T0 a0, T1 a1, T2 a2, T d, T X1) {


  cs[0] = a2[0] / a1[0];
  ds[0] = d[0]/ a1[0];
  for(int i = 1; i < n-1; i++){

    cs[i] = a2[i] / (a1[i] - a0[i-1]*cs[i-1]);
    ds[i] = (d[i] - a0[i-1] * ds[i-1] ) / (a1[i] - a0[i-1] * cs[i-1]);
  }
      
      
  ds[n-1] = (d[n-1] - a0[n-2] * ds[n-2])/(a1[n-1] -a0[n-2] * cs[n-2]);
  
  //求解并输出到控制台
  X1[n-1] = ds[n-1];
  //cout << x[n-1]<<endl;
  for (int i = n-2; i>=0; i--){
      X1[i] = ds[i] - cs[i] * X1[i+1];
      //cout << fabs(x[i]) << endl;
  }
}


# endif
